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ABSTRACT 

Two-dimensional Boussinesq convection is studied numerically using two different meth- 
ods: a filtered pseudospectral method and a high order accurate ENO scheme. The issue 
whether finite time singularity occurs for initially smooth flows is investigated. The nu- 
merical results suggest that the collapse of the bubble cap reported by Pumir and Siggia is 
unlikely to occur in resolved calculations. 1 he strain rate corresponding to the intensifica 
tion of the density gradient across the front saturates at the bubble cap. We also found that 
the cascade of energy to small scales is dominated by the formation of thin and sharp fronts 
across which density jumps. 
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§1. Introduction 


In this paper, we present the results of a careful and detailed numerical study of the 
small scale structures in two-dimensional Boussinesq convection m the absence of viscous 
effects. In particular, we address the issue of whether a finite time singularity can form from 
a smooth initial data. Recently in their numerical study of the same problem! 141 , Pumir and 
Siggia observed that the cap of a symmetric rising bubble (with smooth density variation) 
collapses in a finite time. In contrast with their observations, our results suggest that such 
collapse cannot occur in a finite time. The strain rate associated with the intensification 
of the density gradient saturates, implying an exponential decay for the thickness of the 
front. This is reminiscent of the situation found in vortex reconnection! 1 ’ 13 * 151 : when two 
vortex tubes are brought together, the axial strain rate saturates and the core of the tubes 
undergoes enormous deformation to avoid reconnection (in the absence of viscosity). Thus 
the inviscid solution manages to escape from forming a finite time singularity. 

There are two main motivations for the study of the two-dimensional Boussinesq con- 
vection. One comes from the potential relevance of this problem to the study of atmospheric 
and oceanographic turbulence, as well as other astrophysical situations where rotation and 
stratification play a dominant role. The second motivation comes from the fact that from 
a computational viewpoint, this is the simplest among the class of incompressible flows 
which exhibit vorticity intensification. In particular, it is an open question whether the 
baroclinic generation of vorticity leads to a finite time singularity. Such singularities, if ex- 
ist, provide an effective mechanism for the cascade of energy to small scales. This scenario 
also provides a convenient basis for various turbulence theories which assume, in one form 
or another, that the (ensemble average of) rate of viscous dissipation of energy remains 
finite in the limit of vanishing viscosity, implying the occurrence of finite time singularities 

in many Euler flows! 121 . 

There is a well-known analogy between the two-dimensional Boussinesq convection 
and the three-dimensional axisymmetric swirling flows: 

( V t + U v r + wv z + y UV = 0 

(1-1) jw, + UUJ r + WUJ Z — r UU — ~ 0 
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Here u = u e r + v e B +w e z is the velocity, u = u z -w r is the azimuthal vorticity. Comparing 
(1.1) with (2.1) - (2.4), we see that the centrifugal force plays a role similar to gravity, and 
the azimuthal circulation plays a role similar to density. Grauer and Sideris! 9 ] were the 
first to seek finite time singularities in this restricted class. Although their numerical result 
still remains inconclusive, it has stimulated a lot of recent work along the same direction, 
including the present study. 

A related problem was studied earlier by Childress! 7 ] where he considered nearly two- 
dimensional Euler flows and derived effective equations governing their slow variation using 
contour averaging methods. Under a change of variables, the effective equation takes a form 
similar to the axisymmetric Euler equations with nonstandard connection between circula- 
tion, radius and the angular velocity component. Childress went further to study numeri- 
cally a simplified version of his effective equations and observed finite time singularities for 
the simplified model. This was interpreted as a signal for the re-three-dimensionalization 
of the original (nearly two-dimensional) flow. 

This paper is organized as follows. In §2. we formulate the problem and present 
some preliminary mathematical remarks. In §3. we describe the numerical methods used 
to study this problem. The numerical results are presented in §4. Some discussions and 
concluding remarks are made in §5. 


§2. Formulation of the Problem: 

The equations describing Boussinesq convection are the following: 


{ Pt + u • Vp = 0 

u, + u • Vu + Vp = - (°\ 

V • u = 0 W 

Here p is the density (this should be the temperature and denoted by 6 or T. But we are 
accustomed to calling it density, and therefore denote it by p) , u = (u, v) is the velocity, 
p is the pressure. We have normalized the gravitational constant to be 1. 
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It is convenient to write (2.1) in the streamfunction - vorticity formulation: 

( p t + u • Vp — 0 

(2.2) < + u • Vcj = - p x 

l — Axp = ui 

Here uj is the vorticity, ip is the stream function: 

(2.3) u = v x - u y u — ip y , v = —ip z . 

Introducing the material derivative Dt = d t + u.V, (2.2) becomes simply 

(2.4) D t p = 0, D t u - - p x . 


It is straightforward to show that if the initial data u(x, y, 0) = u 0 (x, y) and p(x , y, 0) 
= po(x,y) are smooth enough, then the solution to (2.1) exists and remains smooth for a 
short time. It is not clear whether any solution would loss its regularity at a finite time. 
However, following Beale, Kato and Majda-'^ , one can show that if a solution losses its 
regularity, then the density gradient has to blow up. More precisely, we have: 


Theorem: Define the norms: 





2 \ 

ll/(0Hm = 

£ / 

dx ai dy a2 
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|/(-)|oo = max \f(x,y)\. 

(x,y)€R 2 

Assume that for some m > 2, ||u(-,0)|| m + ||p(-,0)|| m is finite, but there exists a T* such 
that ||u(-,T*)|| m + ||p(-,T*)|| m = +oo. Then 

• T* 


(2.5) 

a nd 

(2.6) 


J 0 

f T f 

Jo Jo 


, dt = + 


oo 


ds dt — -f oo. 


Proof: We will only give a sketch of the proof since it follows closely the argument 
of Beale, Kato and Majdat 3 !. We will use C to denote a generic constant. 
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Straightforward energy estimates lead to 


Jt IK, Oil™ < !Vu(-,i)U IK, Oil 2 ™, 
|ll»(-,<)llm < |Vu(-,()U ||u(.,<)|| 2 m + IK, Oil 

The logarithmic Sobolev inequality gives us, for m > 2 


Therefore, if we let 


we have 


|Vu(-,t)|oo < CK-,<)|oo(l + log ||u(-,t)|| m ). 


V(t) = \\p(;t)\\ 2 m + l|u(.,OII», 


V(t ) < C|w(-,i)loo (1 + log y(t))y(t). 
Solving this differential inequality, we get 

' f M ,*)lc 

Jo 


y(t) < e e 

On the other hand, since Dtu) — —p x , we have 

d 


y(o)- 


di 


«(',*)|oo < \Px(-,t)\ c 


Hence 


and 


w 


•)0loo < f |/0 X (-, s)|oo ds + C 

Jo 

C f f \px (' i Oloo ds 

p J o J 0 


y(t) < e e - ^ y( 0). 

Since y(0) is finite, we conclude that if y(T*) = +oo, then (2.5) and (2.6) hold. 


§3. The Numerical Methods: 

Computing singular or nearly singular solutions is a difficult task, particularly so in 
the context of incompressible flows. To obtain maximum information one has to push the 
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numerical method to the point where the flow is only marginally resolved. In this situation 
any numerical method is likely to exhibit its own artifact. Such numerical artifacts may 
not go away under simple mesh refinement checks, and they do not necessarily manifest 
themselves in scales comparable to the grid size. Therefore, there is a real danger of being 
mislead by a particular numerical result. To avoid this, we have solved (2.1) using two 
different numerical methods: a Fourier - collocation method and an ENO scheme. When 
the solution develops large gradients, the Fourier - collocation method usually exaggerates 
the situation, whereas the ENO scheme tends to smear out the large gradients. By com- 
paring the numerical results using the two different methods, we can better judge what 
phenomena is likely to be physical. Below we will describe separately the two methods. 
3.1 Spectral Method with Smoothing 

This is the standard Fourier-collocation method^ with smoothing or de-aliasing. 
Roughly speaking, the differentiation operator is approximated in the Fourier space, while 
the nonlinear operations such as multiplications are done in the physical space. We used 
the intrinsic Cray FFT routines which considerably enhanced the performance of the code. 

Since the solutions have large gradients, it is crucial to add filters to the spectral 
method in order that the numerical solutions do not degrade catastrophically if some part 
of the steep gradients are not adequately resolved. A robust way of adding the filters^ 0 ) 
is to replace the Fourier multiplier ikj for the differentiation operator d Xj by 
where 

k mf 

(3- 1 ) <p(k) = e ^ , for |&| < N. 

Here N is the numerical cutoff for the Fourier modes, m, is the order of the filter, and o 
is chosen so that <p{N) = e~ Q = machine accuracy. The machine accuracy on Cray-YMP 
with single precision is roughly 10 14 . Denote by F and F~ ] respectively the forward and 
backward Fourier transform operators, then the numerical derivative is evaluated as 

( 3 - 2 ) D N f = F-\ikip(\k\))Ff. 
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The accuracy of such an approximation scheme depends on the parameter m/. For smooth 
functions /(a:), we have 

(3.3) \\f(x)-D N f(x)\\=0(N- m ') 

We will denote the Fourier- collocation method with rnf - th order filter as SPmf. Unless 
otherwise stated, the results presented in § 4.2 were computed with m f = 10. 


3.2 The ENO Scheme 

The ENO (Essentially Non-Oscillatory) scheme is used for the convection (spatial) part 
of the flow. We use ENO schemes based on point values and numerical fluxes, developed 
by Shu and Osher in [16], [17]; see also [18]. 

To apply ENO approximations, the equation is first written in a conservation form. 
For example the first equation in (2.2) is written as 

pt + (up) x + (vp) y = 0 ( 3 - 4 ) 


The ENO operator is then applied to each of the conservative derivatives in a dimension- 
by-dimension fashion: when approximating (up) x , y is fixed. Unlike the compressible flow, 
the incompressible flow equations are naturally written in characteristic form. Thus no 
expensive characteristic decomposition is needed. Upwindmg can be simply determined 
by the signs of u and v. 

We shall only briefly describe the approximation of a single derivative, say f x . More 
details can be found in [16], [17], and [18]. 


The conservative approximation is of the form: 


(/*);—— (f i+k 


(3.5) 


which, for r-th order ENO scheme, approximates the derivative to r-th order: 




(3.6) 
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where the numerical flux is obtained by interpolating the point values of / on a stencil 
of r + 1 consecutive grid points. The stencil is chosen inductively as follows. For r = 1, 
we choose the stencil to be [j-1, j] or [j, j+1] depending on the sign of u. For r > 1, left or 
right neighboring points are added to the stencil at the previous level v — 1 according to 
the absolute value of the divided differences they each give. In most of our calculations, 
we used r = 3. We will denote this method by EN03. Notice that the scheme is actually 
r + 1-th order in L \ norm. That is, the third order ENO scheme we use is actually fourth 
order accurate in the L\ norm. 

The potential equation in (2.2) is solved with a fourth order central differencing (im- 
plemented in the Fourier space via FFT) plus a fourth order exponential filtering in the 
Fourier space described above. This guarantees fourth order accuracy and in most cases 
avoids instability. We are not sure whether ENO should or can be applied to this potential 
solver, since it is an elliptic equation with possibly a singular right-hand-side. Unlike the 
compressible case, some small oscillations can still be seen in the ENO solution, probably 
due to this potential solver. 

For the temporal discretization, we used Runge-Kutta methods of various order de- 
signed ini 17 ). No major difference between the 3rd, 4th and 5th order methods were found 
in the numerical results. It seems to be a general fact that temporal accuracy is much less 
important than the spatial accuracy. We used the third order version most often since it 
only requires three auxiliary arrays, whereas the fourth order version requires five auxiliary 
arrays. We take initial data that is periodic with period D where D = [0, 2ir] x [0, 27 t], The 
results reported below were computed using CFL equal to 0.5. This is very much within 
the stability region of these methods. 

We point out that our experience strongly favors the use of high order schemes be- 
cause of their small dispersive errors. In a marginally resolved situation, even though the 
filtered spectral method does generate small oscillations near the front where the density 
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experiences a jump, the oscillations are effectively localized near the front. This is because 
the spurious numerical oscillations travel at roughly the right speed which is the speed of 
the front. The higher the order of the filter, the more accurate the propagation speed, the 
more localized the oscillations. Of course the Gibbs phenomena eventually prevails and 
the numerical result becomes noise. 


§4. Numerical Results 

4.1 Formation of the front 


We have numerically integrated (2.1) with a variety of initial data. The initial density 
is usually a perturbation of the uniform state, and the perturbations are localized in each 
period. The early time evolution of such flows is characterized by the formation of a front 
across which density varies sharply. As time evolves, the front gets increasingly sharper 
and the tail starts to roll up. Figures 1-3 present the time development of such events 
for the initial data 


(4.1) 

where 


w(x,y, 0) = 0 

p(x,y, 0) = 50pi(x, y)p 2 (x, y) (1 - />i(x,y)) 


p,(x,y) =( ex »( 1 - + (!/ - *) 2 < * 2 

[ 0 otherwise, 

p 2 (x, y) = { ex P i 1 - ( 1 JKZ- 2 ^ ) ’ if |x — 27r| < 1.95* 

l 0 otherwise 

Figure 1 is the density contour at t = 1.6 computed using SP10 on a 512 2 grid. At this 
time, the flow looks roughly like a rising bubble. Figure 2 presents the same information 
at t = 2.5. By now the outer boundary of the bubble has become a sharp front. We notice 
that as the bubble rises, it leaves behind a long and thin filament of light fluid. This is a 
check on the amount of numerical diffusion present in the scheme. A low order method 
with numerical viscosity (needed to stabilize the front) will destroy the thin filament. In 
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Figure 3a and 3b we present respectively, the contours of density and vorticity at t — 2.9. 
Now the filament of light fluid is thinner and longer, and the trailing edge of the front 
has rolled up. From Figure 3b we see that most of the vorticity resides on the front, the 
rest of the flow domain has very little vorticity. This fact is more drastically displayed in 
Figure 4a and 4b where we plot respectively, the slice of ca and p x at y = n. Since u> and 
—p x have the same sign, from the vorticity equation (2.2) we see that the vorticity peaks 
will be increasingly and monotonically sharper. However, it appears that this mechanism 
of baroclinic vorticity generation only gives rise to an exponential growth of vorticity. We 
will come back to this point later. 

To better appreciate the nearly singular behavior and the computational difficulties 
we face, we will show a few more plots for the profiles of the density and the velocity. 
Figure 5a displays the evolution of the density along the symmetry axis x = 7r at times 
t = 0.5, 1, 1.5, 2 and 2.5. It shows clearly the formation of a front, similar to the formation 
of shock fronts in the solutions of the Burgers equation. At t = 2.5, SP10 on the 512 2 grid is 
not resolving the flow and small numerical oscillations appear on the profile. However the 
numerical oscillations die out on the refined grid 1024 2 . The result of the latter calculation 
is superimposed in Figure 5a on the solution of the 512 2 grid. In Figure 5b and 5c we 
show the slices of p and v respectively, at y = 7r, t — 2.5. From these graphs one might be 
tempted to conclude that jumps have occurred in p and cusps have occurred in v. However, 
a closer look suggests that the cusp-like behavior in v and the discontinuous behavior in 
p go away under refinement. In Figure 5d we compare the numerical results for v on the 
512 2 and 1024 2 grids in the cusp-like regions. It is seen that on the 1024 2 grid the solution 
looks much less singular. Similar phenomena were observed for p. 

These results, particularly Figure 2 and 3a motivate the following question: Is it 
possible that the solution losses its regularity at the same time at the entire front? For 
this to happen there are only two possibilities. One is that while v develops cusps, p stays 
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smooth across the front. This is ruled out by Theorem 1 since it dictates that p and u 
must become singular at the same time. The second possibility is that p develops jumps 
across the entire front. In this case it is necessary that fluids with intermediate densities 
are suddenly swept to the back of the bubble. Clearly this requires infinite velocity whereas 
all of our numerical results give velocities that are very well bounded. 

Before ending this subsection, we present the comparison of the numerical results 
computed using the two different methods described in §3. Figure 6a and 6b display 
respectively, the time history of the maximum and minimum density computed using SP10 
on 256 2 , 512 2 , and 1024 2 grids and EN03 on the 512 2 grid. These quantities should be 
conserved by the exact solutions. As expected, ENO does a much better job in avoiding 
overshoots and undershoots. Figure 6c compares the numerical results for p at y = 7T, t = 
2.5, computed using SP10 and EN03 on 512 2 grid. The result of SP10 undershoots much 
more than the result of EN03, although the latter also contains some small numerical 
oscillations. On the other hand, although not shown here, ENO is usually less accurate 
than the Fourier-collocation method and contains more numerical dissipation. 

4.2 Evolutions of the bubble cap and other aspects of the flow. 

When the bubble rises, lighter fluid has larger acceleration and heavier fluid has smaller 
acceleration. This results in the formation of fronts. However, once the front is formed, 
the pressure gradient becomes important across the front. In fact, as is shown in Figures 
7a and 7b, the formation of increasingly larger pressure gradients may reverse the initial 
picture. In Figures 7a and 7b, we plot the slice of — p y , p , and the acceleration — (p + p y ) 
along the symmetry axis x = tt at two different times t = 0.5 and t = 2. At t = 0.5, the 
acceleration — {p + Py) decreases across the front, whereas at t = 2 it increases across the 
front. (The density p always decreases across the front, see Figure 5a) This implies that 
the velocity difference between the fluid particles at the tip and the back of the front will 
not increase further after t = 2, assuming that the picture remains valid. 
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From another viewpoint, the fluid on the symmetry axis satisfies 

( Pt + vp y = 0 

( 4 . 2 ) < u = 0 

l Vt + Wy =-(/> + Py). 

Let i ^ = p y , p = v y . Then we have 

( 4 . 3 ) £t + v£ y = —t)£. 

The growth of £ is controlled by the strain rate y. In Figure 8a we exhibit the time history 
of the maximum of \p\. During the time interval displayed, the maximum of |r/| always 
occurs at the bubble cap. The results of several different calculations agree remarkably 
well; the strain rate ?/ saturates at about t = 1.7. Consequently, the maximum of |£| can 
at most grow exponentially. This is indeed so, as can be seen from Figure 8b. Other 
quantities, such as u and p x , also grow at an exponential rate. 

We have also done calculations with other sets of data. Here we will briefly report the 
results of a few cases. 

Figures 9a - 9c. display the results for initial data 

u(x,y,0) = 0 

{ . (i-ir) 2 +(»-ir ) 2 

(2tt — x)e 2 5 2 -(*-») 2 -( i /-’ 0 2 

if (x - t r) 2 + (y - tt) 2 < 2.5 2 
0, otherwise. 

The contour curves of p at t = 4.5 and t = 5 are shown respectively in Figures 9a and 
9b, whereas Figure 9c shows the time history of the strain rate \r/\ at the bubble cap. We 
observe that the strain rate begins to saturate at t ~ 4.25 and then increases quickly at 
t Rs 4.8. This corresponds to the time when a mushroom forms out of the bubble cap, and 
the front changes its global configuration. We see that eventually the strain rate saturates 
again at about t = 5.25. We expect that the strain rate at the cap will not increase unless 
the flow changes its geometry at the cap. 
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Figures 10a - 10c present the numerical results for a set of more complicated initial 
data. The data is designed so that at later times the bubble front will develop two bumps 
instead of one as in Figure 9b. Again the increase of the strain rate at the cap is accom- 
panied by a change of the front configuration: after t = 3.5, the front is no longer convex 
(as it was before t = 3.5). We were not able to carry this computation further to see the 
eventual saturation, since the side arms of the bubble become so thin that they cannot be 
resolved on a 1024 2 grid. It is not clear whether a pinch actually occurs at the side arms. 


§5. Discussions and Conclusions 

Our numerical results strongly suggest that the collapse of the bubble cap observed by 
Pumir and Siggia^ 14 ! is unlikely to occur in a resolved calculation. The increase of the strain 
rate is always associated with a global change of the flow character. The contribution of the 
local strain tends to saturate. The global character of the flow is important. Consequently 
any numerical procedure which employs some sort of local approximation is open to serious 
suspicion. This includes the ones presented int 14 i. However, we did not address the issue 
whether other types of singularities can occur. 

Let us reiterate the argument presented in §4.2. Along the symmetry axis x = i r, we 

have 


(5.1) 


f Dtt = -r,£ 

1 D tV = ~{py + Pyy) ~ P 2 


where £ = p y , rj = v y are both negative at the front. From (2.1) we also have 


(5.2) -p y - Ap = u\ + V 2 + 2 UyV X . 

On the symmetry axis, this reduces to 


(5.3) 


Py Pyy — Pxx T 2 rj . 


12 


Therefore, we have 


(5.4) 


Dtt = - Vi 
Dtn = Pxx + v 2 - 


As long as the front remains relatively smooth, p xx is negligible compared to 7/ 2 . Conse- 
quently £ and // will stay bounded. In fact since near the cap p y < 0, p yy < 0, we obtain 
a necessary condition for collapse at the cap: 


(5.5) -2 r/ 2 < p xx < -T] 2 . 

In particular, the front itself has to develop infinite curvature. This is consistent with the 
argument presented in section 5 oft 14 !. 

The key difference between the computations reported ini 14 ) and the computations 
reported here, is that ini 14 ) a high frequency instability (with a scale that is comparable to 
the thickness of the front) is observed at the cap, whereas in all our resolved calculations 
we have not seen the manifestation of such an instability. Since the front is extremely thin, 
near the cap the situation resembles that of the Rayleigh- Taylor instability. Therefore a 
linear stability theory will indeed predict that the front is unstable to short wavelength 
perturbations. However, the relevant question here is not whether perturbations can grow, 
but whether short wavelength perturbations can exist in a dynamic context. In other 
words, the issue is whether the exact solutions of (2.1) can supply perturbations with 
wavelength on the scale of the thickness of the front. 

We want to emphasize that the question we are asking is a rather academic one. In 
actual physical experiments, various external noises provide ample sources of such pertur- 
bations. The solution we are seeking numerically may hardly occur in experiments. (This 
is yet another argument for the importance of scientific computing: It has the potential to 
carry out much more controlled experiments.) 

Before discussing the question raised above, let us recall two problems that bear 
some similarity with the case we are studying. The first is the classical Kelvin- Helmholtz 
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instability of a vortex sheet. The braid that connects two consecutive rolls in a Kelvin- 
Helmholtz roll-up is still subject to Kelvin-Helmholtz instability, although the local vortex 
sheet strength is decreasing in time and the braid is under straining. This kind of instability 
manifests itself in numerical computations in such a way that the round-off error causes 
small spurious roll-ups along the braids. This phenomenon is studied very carefully by 
Krasnyt 1 d and is cured by introducing a nonlinear filter that keeps the numerical solution 
effectively analytic. 


The second problem is the behavior of an expanding bubble in a Hele-Shaw cell with 
very small surface tension. The presence of the surface tension makes the interface between 
two viscous fluids linearly stable. Nevertheless, the interface is nonlinearly unstable and 
prefers to develop the characteristic viscous fingers and the tip-splitting of these fingers 
observed in experiments. For example the amount of perturbation needed to upset the 
linear stability of a viscous finger and cause it to split is exponentially small in terms of 
the surface tension^ . In the numerical studies of Dai and Shelley^ , the lower precision 
calculations give rise to finger-like behavior. One might argue that such fingers are physical, 
on the basis that they are observed in experiments. But they are not faithful solutions 
of the original dynamic problem. Indeed as was shown by Dai and Shelley, the initial 
instability go away in the high precision calculations (128 bits). 


Returning to the problem we are studying, the linear stability of the cap of a rising 
bubble is investigated by Batchelor et. alJ 2,14 h In particular, Pumir and Siggia show that 
perturbations with wavelength 8 (the thickness of the front) can amplify by a factor of 
(where -y c is the radius of curvature at the cap) if 


(5.6) 


7c 1 

6 28s 2 


2 ( 7-1 

< e ^277^ 


> 1/2 


-i) 


Here s is the strain rate at the cap. We remark that an infinitesimal disturbance eventually 
decays since the local strain increases its wavelength exponentially. 
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Although in some cases we have carried out the computation well beyond the limit 
set by (5.6), we have not observed the kind of instability at the cap reported by Pumir 
and Siggial 14 !. One such instance is displayed in Figure 11 where we have roughly 8 = 

j c = 1,3 = 2. (5.6) becomes 2048 < 1.109 x 10 4 and is very well satisfied. Still no 
instability on the front is observed. 

A similar situation occurs in the classical Rayleigh- Taylor problem. The interface 
between a heavy fluid on top and a light fluid at the bottom is linearly unstable and 
typically evolves to mushroom-like or plume-like structures. The tips and tails of the 
structure are still subject to linear instability. But no instability at the tips or tails has 
been reported, even though the problem has been extensively studied^ 10 ’ 19 !. 

We remark that for the collapse to occur, a cascade of instabilities has to develop at 
the cap with increasingly smaller temporal and spatial scales. We have argued that such 
an instability is unlikely to arise for the true exact solutions of (2.1). However, since the 
front is very sensitive to numerical errors, any uncontrolled truncation can easily trigger 
the instability and lead to erroneous results. 

With respect to the issue of formation of small scale structures, our calculation sug- 
gests that the cascade of energy to small scales is dominated by the formation of thin 
fronts. Roll-ups that are characteristic of two-dimensional Kelvin-Helmholtz instability 
also occur but play a less important role. Occasionally we also observe small rolls formed 
out of big rolls, evidence of what is suggested in Richardson’s well-known quote. But these 
rolls carry a very small amount of energy. The scenario of a cascade through a sequence 
of roll-ups on increasingly smaller scales, if at all possible, has negligible contribution to 
the overall cascade mechanism. 

The evolution of the energy spectrum for the first set of data is plotted in Figure 
12. There is an apparent similarity to the spectrum observed by Brachet et. al. for the 
Taylor-Green vortex^. The near singular behavior creates an algebraically-decaying part 
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in the bulk of the spectrum. However the spectrum does not seem to settle down to any 
stationary state. 

To conclude, let us emphasize that the ultimate answer to the question addressed here 
must be analytical. Numerical studies can at best provide partial evidence. For incom- 
pressible flows, it has proven very difficult to obtain analytical results on the regularity of 
the solutions. As a result, considerable efforts were devoted to the numerical approach. 
However, computing singular or nearly singular incompressible flows is also an extremely 
difficult task. Not surprisingly, the subject of seeking singularities numerically for incom- 
pressible flows is full of controversy. The contradicting conclusions reached in the present 
paper and the work of Pumir and Siggia^ 14 ! certainly adds one more to the list of many 
unsettled issues. It is our intention that these numerical work will generate sufficient inter- 
est for further analytical study and the evidence presented will contribute to the ultimate 
resolution of the issue raised. 
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Figure 1. Density contour at t=1.6 for initial data (4.1). 
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Figure 2. Density contour at t=2.5 for initial data (4.1). 
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Figure 3a. Density contour at t=2.9 for initial data (4.1). 
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Figure 3b. Vorticity contour at t—2.9 for initial data (4.1). 
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Figure 5a. Evolution of density along the symmetry axis x = 7r at t - 0.5, 1, 1.5,2, 
and 2.5, computed using SP10 on a 512 2 grid. The result at t = 2.5 computed on a 1024 2 
grid is also superimposed. Notice that the small numerical oscillations disappear on the 

refined grid. 



x 

Figure 5b. Slice of p at y — 7r, t = 2.5. 
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Figure 6a. Time history of the maximum density from different computations. 
From top to bottom: SP10 on the 256 2 grid, SP10 on the 512 2 grid, SP10 on the 1024 2 
grid and EN03 on the 512 2 grid. 



Figure 6b. Time history of the minimum density from different computations. 
From bottom to top: SP10 on the 256 2 grid, SP10 on the 512 2 grid, SP10 on the 1024 2 
grid and EN03 on the 512 2 grid. 
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Figure 8a. Time history of the maximum of |^| at the front (which is also the 
strain rate v y at the cap). From top to bottom: SP10 on 1024 2 grid, SP10 on 512 2 grid, 
and EN03 on 512 2 grid. Strain rate saturates at t = 1.7. 



t 

Figure 8b. Time history of the maximum of |£| . From top to bottom: SP10 on 
1024 2 grid, SP10 on 512 2 grid, and SP10 on 256 2 grid. 





Figure 9a. Density contour at t = 4.5 for the initial data (4.4). 
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Figure 9b. Density 


>ur at t = 4.5 for the initial data (4.4). 
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Figure 10a. Density contour at t = 3.5 for another set of initial data. 
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Figure 10b. Density contour at t = 5 for the same data. 
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Figure 10c. Time history of the strain rate at the cap. 
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Figure 11. Contour plot of density in another set of calculation. 
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